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LOW DEGREE RESONANT GEOPOTENTIAL 
COEFFICIENTS FROM EIGHT 
24-HOUR SATELLITES 

C. A. Wagner 

ABSTRACT 

The most extensive tracking data base on 24 hour satellites has been analyzed 
for resonant geopotential effects in the period 1963-1969. The observations used 
consisted of 279 sets of mean Kepler elements in 21 free drift arcs of eight sat- 
ellites, each arc between 7 weeks and 3-1/2 years in length. 

Current geopotential solutions from the Smithsonian Astrophysical Observa- 
tory (SAO) which do not use the 24-hour data are shown to recover most of 
these observations almost as well as in solutions using this data. However, 
since the longest arcs are poorly recovered with the "best" Smithsonian field, 
an improved resonant geopotential is desirable. Therefore, a new solution is 
derived which combines the best available low degree Smithsonian coefficients 
with the 24 hour data. This solution is virtually indistinguishable from the "best 
fitting" solution from only 24 hour data. A comparison of this combined solution 
with the other current solutions shows that: 1. the east-west equilibrium points 
for geostationary satellites are now known accurately to better than 1/2 degree, 
and 2. the resonant acceleration on a 24 hour satellite can now be predicted ac- 
curately to better than 10' 6 radians/day 2 . 
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The computational efficiency of this resonant solution should be noted. 
The solution essentially uses only 279 pieces of data compared to over 50,000 
in the SAO geopotentials: the resonant coefficients themselves agreeing re- 
markably well. The total processing time for the resonant solution (involving 
five iterations) was less than 4 minutes on an IBM 360/91 computer. 



LOW DEGREE RESONANT GEOPOTENTIAL 


COEFFICIENTS FROM EIGHT 
24-HOUR SATELLITES 


INTRODUCTION 

Since 1963 more than forty 24-hour satellites have been placed in orbit by 
the National Aeronautics and Space Administration and the Department of De- 
fense.* Of these, only about twenty have been close enough to exact commen- 
surability with the earth's rotation to be especially useful for the purposes of 
dynamic satellite geodesy; determining the geopotential from measuring the 
perturbations of satellite trajectories. The tracking record of the eight satellites 
discussed here represents well all the geographic longitudes covered by all the 
nearly commensurate objects for which data have been available (i.e. unclassi- 
fied) from 1963 to 1969. 

The theory of the greatly amplified long term perturbations of these synchron- 
ous, or resonant, orbits has been much discussed in the literature. 3 The analysis 
presented here is based on the rapid numerical evaluation of perturbations to 
mean Kepler elements over long tracking arcs 3 (7 weeks to 3-1/2 years). The 
observations used nre themselves mean Kepler elements determined from "raw 
tracking data" over short (generally 1-2 week) arcs by outside orbit programs of 
varying degrees of sophistication. The ultimate "raw" tracking data in these arcs 
are: range and range rate and minitrack (radio interferometric angles for the 
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Synchronous Communication Satellites (SYNCOM’s); radar range, azimuth and 
altitude for the International Telecommunication Satellites (INTELSATS); and 
range and range rate for the Applications Technology Satellites (ATS's). 

The use of Kepler element data from deeply resonant satellites to aid in 

the definition of the longitude dependent geopotential has also been extensively 

discussed in the literature.^ - ® Previous results for 24 hour satellites have 
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depended on much less data than is available now. 

The methods of analysis used previously have also been somewhat limited 
so that the full resonant information content in the Kepler element evolution 
has not been used. For example Allan and Piggott^ in 1966, used the most 
sensitive longitude information (the equator crossings) directly, fitting this data 
to a numerically integrated trajectory. But the sensitive information in the 
semimajor axis was ignored in their results. On the other hand Wagner*®, in 
1968, used both satellite position and orbit energy but indirectly, to derive ac- 
celerations which were fitted to a purely analytic model. However, by far the 
greatest deficiency in all previous results has been due to the sparse longitude 
coverage of the data. This deficiency has been greatly aleviated by the addi- 
tion of ATS and INTELSAT data filling large gaps in the previous record. Con- 
tinued tracking of SYNCOM 2 and 3 has also made valuable extensions to the 
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ANALYSIS 


The analysis here is an application to these orbits of the geopotential 
recovery method described in Wagner (1969). 3 Two points of departure should 
be mentioned. 

1. In this analysis, the partial derivatives of Kepler elements with respect 
to the initial arc parameters and geopotential coefficients are calculated by 
numerical integration of variation equations derived from the Lagrange planetary 
equations. 

The Lagrange equations (for the evolution of a perturbed trajectory) are 
abbreviated by: 


e . 


f i (<V c k ) , i, j = 1, 2, 


k = 1, 2, 


where 

n = number of gravity constants solved for. 
e . = Kepler element. 

Then the "variation equations" for the c k constants ai'e found by taking the partial 
derivatives of these Lagrange equations with respect to the c k , s (using the chain 
rule) and interchanging total with partial differentiation: 



i - 1. 2, ••• 6 - 


3 


% 


\ 


The bracketed quantities are the required "partials" necessary for the first 
order differential correction of the c k 's with respect to observations of the e. . 
The Kepler elements analyzed in this fashion are the classical set: semi- 
major axis, eccentricity, inclination, argument of pericenter, right ascension 
of the ascending node and mean anomaly. 

The so-called "state transition" equations (for the 6 x 6 = 36 partials of the 
set of e with respect to the set of e Q , the initial Kepler elements in the arc) are 
easier to write down because there is no functional dependence here of the gravity 
coefficients with respect to the initial Kepler elements. Strictly speaking, of 

ft 

course, there is such a dependence because the data used are mean elements. 

But the observations here are given as fixed (i.e. not corrected for gravity model 
changes inside or outside this long arc orbit determination process). Thus there 
is no formal dependence. In order to introduce such a dependence, it would be 
necessary to go back to the original observations and use the corrected gravity 
field to produce new mean elements. The changes so produced would be entirely 
insignificant compared to the changes in the long arc reflecting the resonance. 

The 36 "state transition" equations, numerically integrated, are: 
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The initial conditions for the complete set of variation equations are simply: 


0. i ? j 
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From the simultaneous integration of these "variation equations" along with 
the trajectory (Lagrange) equations, the observation partials and residuals are 
found and normal equations accumulated for a standard weighted least squares 
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process. For example, the observation equations formed for the first order 
differential correction process, using these partials, are: 


. (observe 


w 

d) - e, (computed)] = £] [^] •'-‘Vi * £ [W] 


An estimate of the accuracy of each observation quantity is used as the weight 
for the corresponding equation of condition (or observation equation) (see Chapter 
5 in Kaula's text 2 ). The programming system which does the differential cor- 
rection process is called ROAD (for Rapid Orbit Analysis and Determination). 11 
ROAD can also accept a-priori information in the form of initial diagonal 
covariance matrices for the reference gravity field and the starting arc elements. 


The inverse of these initial covariance matrices forms the initial estimate 


of the normal matrix. The initial normal vector of residuals, the right hand 
side of the normal equations, is zero. It is recalculated at each subsequent 
iteration as the product of the a-priori normal matrix and the difference between 
the new and a-priori solution vectors. 

2. In the present 24 hour solution, the only data actually fit -to was the mean 
longitude of the ascending equator crossing (also called the mean longitude) which 
is physically the most meaningful and sensitive resonant orbit parameter. 12 This 
quantity, for 24 hour satellites, is the Greenwich hour angle subtracted from the 
sum of the Kepler elements mean anomaly, argument of perigee and right ascen- 
sion of the ascending node. The semimajor axes, though available were not used 
in the analysis because their definitions vary from arc to arc. For example, in 
some arcs the luni-solar perturbations in the semimajor axis have been removed 
analytically or by smoothing (where mean Kepler elements have been originally 
reported,. In other arcs these perturbations have been left in (where osculating 
elements have been originally reported). The present ROAD preprocessor converts 
osculating Kepler to mean Kepler elements, but only accounting for the effects of 
geopotential zonal harmonics using Brouwer's theory. 12 This conversion is ac- 
complished by subtracting from the given osculating elements, the short period 
perturbations (due to the zonal harmonics) calculated from Brouwer 5 s theory. 

It is difficult to find a definition of mean semimajor axis completely con- 
sistent with the ROAD integration even in a single arc. The real problem is in 
the long term evolution of the mean anomaly which is extremely sensitive both 
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to small semimajor axis changes and to the resonance effects that are to be 
determined. The most convenient solution to this problem would appear to be 
to allow the ROAD program to correct automatically for a small constant shift 
in the semimajor axis observations for each long arc. In this way the semimajor 
axis data could be utilized simultaneously with the longitude data without the need 
for manually correcting the observations. The statistics (or weight) of the 
present solution should be improved considerably when the semimajor axis data 
are used. However, major improvement in the results are only expected when 
additional longitude gaps are filled. 

DATA AND PREPROCESSING 

The actual data used (for fitting purposes) in this solution was the mean 
longitude of the ascending equator crossing (lambda, \) as described previously. 
However, in Table 1 1 list the full set of mean Kepler elements available for the 
21 arcs. The Table 1 column headings are: Time is the epoch in modhied 
Julian days (Julian days - 2400000.5), A is the orbits semimajor axis in earth 
radii, E is the orbits eccentricity, Incl is the orbits inclination in degrees, 

Omega is the orbit's argument of perigee in degrees, Node U ti*e orbit's right 
ascension of the ascending node in degrees, Lambda is the orbit's mean geo- 
graphic longitude in degrees. The reference and solution trajectories duplicated 
these elements quite faithfully although they were not specifically designed to do 
so. In future solutions specific use of the other Kepler element data in Table 1 
will be made, to improve this solution. 
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The original orbit elements given from short arc solutions varied from 
satellite to satellite. The source of this original data and it's subsequent 
processing to the mean elements in Table 1, is given below. 

SYNCOM 2 Arc 1 

The original data were mean elements derived from Goddard range and 
range rate and minitrack data, as determined (at NASA -Goddard) from an 
analytic program using Brouwer's orbit theory.13,14 These "Brouwer" mean 
elements were used directly by ROAD in this solution. 

SYNCOM 2 Arcs - 2, 4, 5 and 8; SYNCOM 3 Arc 6 

The same kind of information applies to this arc as that above except that 
essentially no minitrack information is contained in the mean elements. 

SYNCOM 3 Arc 7 

This arc contains the same kind of elements as those above. However 
there was an orbit maneuver 2 weeks after the first set of elements which re- 
duced the inclination by about 1° but did not substantially change the orbit energy 
or longitude elements sensitive to the resonance effect. Nevertheless there has 
been some difficulty getting good individual arc solutions with this data. In the 
future the first two data sets will probably be dropped when many arcs are 
processed simultaneously. 
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SYNCOM 2 Arc DoD, SYNCOM 3 Arcs 11, 13 and 14 


The original data for these arcs were osculating polar coordinates of the 
spacecraft: Right ascension, declination, distance from the earth's center; 
azimuth, elevation and speed of the velocity. These elements were reported 
by the Air Force Systems Command, Sunnyvale California, from range and 
range rate data. The numerical orbit theory used to derive these elements in- 
cluded sun and moon gravity effects and radiation pressure as well as a fairly 
accurate set of fixed geopotential constants including resonant ones. However 
in no case did the data span for each set of elements extend beyond 3 weeks. 
Thus the orbit programs resonant constants served mainly to define more ac- 
curate individual elements for the spacecraft, for prediction purposes. 

These osculating polar elements were converted to osculating Kepler ele- 
ments^® and then to mean Kepler elements using Brouwer's theory.13 Thus 
they do not have the short (less than 1 day) period gravity effects of the sun and 
moon removed from them. The most serious result of this deficiency is the 
introduction of a significant constant bias in the semimajor axis which makes 
that data incompatible with the evolution of the longitude data when both are 
used in ROAD . 

ATS 1 and ATS 3 Arc 1 

The original tracking data in these arcs were Goddard range and range rate 
observations. The preprocessing consisted of fitting data in these separate long 
arcs to a trajectory computed numerically by the GEOSTAR program. This 
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program found individual best fitting values of the harmonics C 22 and S J2 ac- 
cording to the data in these arcs. (The program also calculated a full set of 
fixed geopotential as well as the usual luni-solar effects.) 12 sets of osculating 
Kepler elements from the GEOSTAR solution trajectories within the data periods, 
were chosen as representative of the evolution for each long arc. The osculating 
data was then converted to mean data through Brouwers theory. Smoothing 
tests of this orbit data were made with the ROAD (mean element) program. 

They revealed that the short period perturbations not removed by the Brouwer 
theory (due to the sun and moon for example) provided "pseudo noise" to the 
GEOSTAR elements of about the same magnitude (0.002° in longitude) as that 
indicated by the original GEOSTAR solution residuals. Thus, in spite of the 
fact that the GEOSTAR solution elements for these arcs can not be proven to 
be strictly uncorrelated, they appeared to be so because of the short period 
perturbations left in them. By this device I was able to produce a set of "noisy" 
elements representative of the excellent range and range rate data over these long 
arcs. It would have been difficult or impossible to obtain equally good "short 
arc" elements (correcting for the resonance effect) in strictly independent data 
spans. There were few periods in ATS tracking when the simultaneous coverage 
from 2 or more ground stations was good enough to obtain such independent orbits. 

ATS 3, Arcs 2 and 3; and ATS 5 

These arcs were treated in the same way as the ATS arcs above except for 
the introduction of artificial noise into the GEOSTAR solution elements. Here 
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it was found that the "noise" introduced because the short period perturbations 
were not removed by the Brouwer theory was too great and too biased to ade- 
quately represent the original GEOSTAR results for these arcs. Therefore an 
additional smoothing was performed followed by a controlled "noising" up of 
. the data, .he procedure was as follows: ROAD (mean element) trajectories 
were fitted to the representative osculating elements in the GEOSTAR solution 
trajectories (using the resonant coefficients calculated in the GEOSTAR long 
arc solutions). Longitudes of the equator crossings from these mean element 
trajectories were then put through a random noise generator and noised up to 
a ± 0.002° (lcr) level (representative of the excellent GEOSTAR results for 
these arcs). The weight actually assigned (in the multi-arc ROAD solutions) 
for all the ATS longitude data was 0.005°, to compensate further for any cor- 
relation remaining from the original GEOSTAR solutions. 
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Early Bird, Arcs 1 and 2 

The original tracking data in these arcs consists of radar ranges, elevations 
and azimuths from the Com munication Satellite Corporations (COMSAT) facility 
in Andover Maine. The large amount of these data were reduced by COMSAT 
to a small set of independent mean Kepler elements with the use of a trajectory 
computation program based on Brouwers 1959 theory but including sun and moon 
gravity effects as well. This mean Kepler data was used directly (without further 
processing) by ROAD in the solutions reported here. 
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I NT F SAT 2F3, Arcs 1 and 2, and INTELSAT 2F4 


For these arcs, COMSAT employed a number of transportable antennas to 
provide Pacific ocean coverage (for INTELSAT 2F4) and assist the Andover 
radar for tracking over the Atlantic ocean. A "Cowell" type osculating element 
orbit determination program (with a limited low degree geopotential field) was 
used to fit the radar data for these arcs. The resulting independent short arc 
"best fitting" osculating Kepler elements reported by COMSAT were converted 
to mean elements by Brouwer's 1959 theory. Once again the short period sun 
and moon effects are not accounted for in this conversion. 

NEW 24-HOUR SATELLITE GEOPOTENTIAL SOLUTIONS 

Many trajectories were determined by the ROAD program for the individual 
arcs in Table 1. These preliminary solutions, including the determination of a 
few critical resonant terms to absorb the residual acceleration in each long arc, 
established the quality or weight of the longitude (\) data. The standard deviation 
of this data type in these individual arc solutions, rounded up to the nearest 0.005° 
to be as conservative as possible, served as the arc weight when the multi-arc 
solutions were made. These conservative standard deviations and other summary 
data from the multi -arc solutions are shown in Table 2. It should be noted in 
Table 2 that many of the arc root mean square (RMS) residuals in the multi-arc 
solutions are indeed smaller than the estimate of the standard deviation from 
the preliminary best fitting trajectories. This is not surprising considering the 
rounding -up of the estimate as well as the greater number of degrees of freedom 
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lost in determining individual best fitting resonance constants for each arc with 
limited data. Indeed, as Table 3 shows the weighted RMS deviation of an ob- 
servation in the best multi-arc solutions is actually less than 1. 

Table 2 presents the root mean square residuals (observed minus com- 
puted values of the mean longitude) in the 21 arcs from 3 of the 4 multi -arc 
solutions whose geopotential coefficients and overall statistics are found in 
Table 3. The multi-arc processing with ROAD had three objectives; illustrated 
by the 3 solutions in Table 2: 

1. To see how well the best nonresonant satellite geopotentials could recover 
the longitude data in fixed field solutions for only the 6 initial Kepler elements 
for each arc. The best such solution was with the 1969 Smithsonian Astrophysical 
Observatory's ( SAP ) COSPAR field. The RMS residuals in the longitude (A.) 
for the arcs in this solution are in the last column of Table 2. The COSPAR field 
coefficients used and the overall weighted RMS residual of this solution is shown in 
Table 3. Except for the especially long arcs SYNCOM 2 DoD and SYNCOM 3-14, the 
COSPAR solution compares satisfactorily with the adjusted coefficient solutions. In 
a number of other arcs however (ATS 5 and INTELSAT 2F4 for example) this best 
nonresonant fixed field solution showed significantly degraded residuals (reveal- 
ing clear systematic errors arising from the geopotential). On the other ham. in 
the excellent ATS 1 and 3 arcs the COSPAR solutions were fully competitive with 
the adjusted solutions (perhaps because these arcs were relatively short). It is 
remarkable that the otherwise excellent SAO 1970 (standard earth II) field,*** 
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which is very close to the COSPAR field in the resonant coefficients sensitive for 


24-hour satellites, gives a solution through this data considerably degraded from 
the COSPAR. (See Table 3.) Yet the accelerations on a geostationary satellite 
calculated from the SAO 1970 field are well within 0.1 x 10" 5 radians/day 2 , or 
almost indistinguishable from those of the preferred adjusted field (see Figure 1). 
One may fairly conclude from these fixed field solutions that only minor adjust- 
ments to them are necessary for use with 24 hour satellites. However one 
should also conclude that even such minor adjustment will produce significantly 
better results with very long 24 hour trajectories (3 months to years in length). 

For this reason I had as my next objective: 

2. To see how much change would be necessary in the latest SAO nonresonant 
fields to recover the longitude data in the 21 arcs to as high an accuracy as pos- 
sible (i.e. in a weighted least squares sense). The first and best such test of the 
24 hour data was made by solving (in ROAD) without constraint for the 3 dominating 
harmonics (2,2), (3, 3) and (3, 1). The arc results of this solution are in column 8 
of Table 2. The field coefficients and the overall statistic of the solution (which 
included fixed SAO 1970 values of (4,2) and (4,4)) are in Table 3. The covariance 
statistics of this solution are displayed in Table 4. The Table 4 statistics can be 
said to be the minimum data necessary to rigorously combine the 24 hour track- 
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ing results with other gravity data for improved geodetic solutions. 

The weak value for (3, 1) and the relatively large number of high correla- 
tions in this unconstrained solution suggested a third objective: 
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3. To see how much improvement in the correlations and the (3, 1) harmonic 
would be possible with a constrained solution through the 24-hour data using the 
1970 SAO field as a-priori information. For this purpose I allowed the following 
lc variations in the resonant SAO field coefficients (unnormalized) through (4, 4) 
(based on the differences among all the 1969 and 1970 SAO geopotentials): (2,2) 
and (3,1), 15 x 10' 9 ; (3,3) and (4,2), 10 x 10' 9 and (4,4), 10 x 10' 10 . With this 
constraint, the weighted least squares multi-arc ROAD solution for all these 
resonant coefficients, using the longitude data, gave the field having the statistics 
shown in Tables 2 (column 9), 3 and 5. Evidently this is a considerably smoother 
and better solution than the unconstrained one as revealed by the lower correla- 
tion coefficients, and the more representative value of (3, 1). In addition, the 
overall arc residuals (Tables 2 and 3) are only slightly degraded over the 
unconstrained solution. 

In conclusion, this constrained solution should be preferred for work with 
24 hour satellites. In Figure 1 is plotted the east-west longitude acceleration on a 
geostationary (E = 0, 1 = 0°, Period = 24 hrs.) satellite as calculated from this 
constrained field and the resonant part of the 1970 SAO geopotential. The formula 
used to compute this acceleration is found in Wagner (1968). Also in Figure 1 
is a histogram giving the number of 24 hour satellite observations used in the 
present solution per 10° of longitude (total = 279) compared to that used in a 
1966 solution by Allan and Piggott. 7 This histogram shows the significant im- 
provement in the data base for 24 hour satellite results in the past 4 years. The 
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biggest gap in data is now between 0° and 60° which will soon be partly filled 
by the tracking of Britain's skynet satellites launched in 1970. 

Table 6 gives a comparison of east-west equilibrium points found from the 
longitude acceleration zero's on the geostationary satellite. The NWL-8D (classi- 
fied) field was determined (by the Naval Weapons Laboratory) from doppler meas- 
urements on the Navy's "transit" (high inclination, low altitude) navigation satel- 
lites. It is recalled that the SAP COSPAR (1969) field (which is used here in the 
constrained solution) resulted from only satellite optical data while the 1970 SAO 
field is relatively independent from it having significant contributions from 
surface gravity and deep space probe measurements. It is clear from this table 
that the equilibrium points can now be said to be known to better than 0.5°. 
Similarly it can be seen from Figure 1 that the geostationary satellite's accelera- 
tions are known to better than 0.1 x 10" 5 rad/day 2 . 

Finally, in Figure 2 are displayed the observations and ROAD computed 
values of the mean longitude (lambda or \ in degrees) in the most revealing arcs 
of the constrained resonant solution. (Only those arcs whose "raw" mean longi- 
tude histories show strong accelerations are displayed. The arcs not shown 
reveal equally significant accelerations when the data plotted is the actual mean 
longitude minus the longitude in a trajectory drifting at the steady average rate 
for the arc.) If no significant resonant perturbations acted on these satellites, 
their longitude histories would be linear. 
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CONCLUSIONS 


Tests of long arcs of 24-hour satellite tracking data (mean elements) in 
the period 1963-1969 reveal: 

1. For most of the arcs less than a year in length, the use of recent 
Smithsonian gravity fields (derived without the use of 24-hour data) can re- 
cover the data satisfactorily. 

2. For the arcs over one year long, the Smithsonian fields cannot adequately 
recover the data due to small deficiencies in low degree resonant geopotential 
terms. 

3. A combined Smithsonian-24 hour data field can recover all the 24 hour 
data with satisfactory precision and is recommended for use with all synchronous 
satellites at any longitude or inclination. 

4. This recommended field should predict the long term geopotential ac- 
celeration on a 24 hour satellite to better than 10" 6 radians/day 2 , and the 
equilibrium points for a geostationary satellite to better than 0.5°. 
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MISSION & TRAJECTORY ANALYSIS OIVISION 

BRANCH DATE 8/70 

BY C.A Wagner PLOT NO. 2178 




RESONANT GEOPOTENTIAL COEFFICIENTS FOR 24 HOUR SATELLITES 
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STRAINED) (1D70) 




RRELATION COEFFICIENTS AND STANDARD DEVIATIONS IN 
UNSTRAINED 24 HOUR SATELLITE GE- 'POTENTIAL SOLUTKJN 




TABLE 6 

EAST -WEST EQUILIBRIUM POINTS FOR 
GEOSTATIONARY SATELLITES 

Gravity Field 

Indian Ocean 
Stable Pt. 

West Pacific 
Unstable Pt. 

East Pacific 
Stable Pt. 

Atlantic 
Stable Pt 

WAGNER-SAO 
COSPAR (CON- 
STRAINED) (1970) 

75.0° 

161.8° 

-105.5° 

-11.5° 

SAO STANDARD 
EARTH n (1970) 

i 

75.1 

i 

162.1 

-105.1 

-11.2 


NWL 8-D 
(1967) 







ACCELERATION USING 24 
HOUR SATELLITES (from 



Figure 1. Longitude Acceleration on a Geostationary Satellite (and Observation Histogram) 
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TIME (MODIFIED JULIAN DATE) 
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(c) SYNCOM 2,2 Data 



TIME (MODIFIED JULIAN DATE) 
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(d) SYNCOM 2, 8 GSFC Data 


Figure 2 (Continued). Observed and Computed Measurements- Lombda 
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(i) SYNCOM 2, 1 Data 


Figure 2 (Continued). Observed and Computed Measurements - Lambda 
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